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We consider two different collective spin systems subjected to strong dissipation - on the same 
scale as interaction strengths and external fields - and show that either continuous or discontinuous 
dissipative quantum phase transitions can occur as the dissipation strength is varied. First, we 
consider a well known model of cooperative resonance fluorescence that can exhibit a second-order 
quantum phase transition, and analyze the entanglement properties near the critical point. Next, we 
examine a dissipative version of the Lipkin-Meshkov-Glick interacting collective spin model, where 
we find that either first- or second-order quantum phase transitions can occur, depending only on the 
ratio of the interaction and external field parameters. We give detailed results and interpretation for 
the steady state entanglement in the vicinity of the critical point, where it reaches a maximum. For 
the first-order transition we find that the semiclassical steady states exhibit a region of bistability. 
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I. INTRODUCTION 

The field of ultracold quantum gases has recently 
made remarkable progress toward the implementation of 
(fully) tunable interacting many-body quantum systems 
p. . Specifically, the degree of control in experiments al- 
lows for a precise variation of system parameters, for ex- 
ample interaction strengths and effective fields, such that 
the systems can be made to undergo transitions between 
different quantum phases [2!. 

Of particular interest are microscopic, interacting 
many-body systems, which have been widely studied in 
the context of closed systems, where quantum phase tran- 
sitions (QPTs) arise due to the competition between fluc- 
tuations originating from different coherent processes in 
a system (e.g. tunneling versus interaction in the Bose- 
Hubbard Model) [3^, '5 . Although individual systems sub- 
jected to dissipation on the same scale as their charac- 
teristic frequency have been extensively studied [5 , the 
effects of dissipation on interacting many-body systems 
are less well-known. Recently a collective spin system 
with weak dissipation was studied in the context of a 
non-equilibrium QPT 0. It was shown that the well 
known second-order phase transition found in the equiv- 
alent closed system persists, with weak dissipation being 
responsible only for minor modifications to the system 
properties. However, in addition, a first-order phase tran- 
sition was shown to occur exclusively due to the presence 
of dissipation, i.e., this phase transition is absent in the 
equivalent closed system case. For both types of tran- 
sition, the spin-spin entanglement was shown to exhibit 
pronounced signatures of the criticality. 

Given these non-trivial results for the case of weak 
dissipation, it is then naturally interesting to consider 
the regime of strong dissipation, i.e., dissipative rates on 
the same scale as the interaction strengths and external 
fields. In this regime marked differences are expected in 
comparison to the closed system case, and, specifically, 
new types of QPTs, driven by the dissipation, are ex- 



pected to emerge [71 H H [TOl [H [Tlll3] . In this work 
we consider two models of open collective spin systems 
where a QPT arises solely due to a competition between 
fluctuations associated with Hamiltonian (coherent) dy- 
namics and with dissipative processes. In addition to 
studying elementary characteristics of the phase tran- 
sitions, we also study entanglement criticality and find 
that pronounced maxima in entanglement measures oc- 
cur at the QPT. A further interesting feature that arises 
in the present work is that, for the second model con- 
sidered, the nature of the phase transition, i.e., whether 
it is continuous or discontinuous, is governed by the ra- 
tio of the spin-spin interaction strength to the effective 
("magnetic") field. This behaviour is in strong contrast 
to the equivalent non-dissipative models, where the char- 
acter of the phase transition is governed by the nature 
of the interaction, i.e., by whether it is "ferromagnetic" 
or "anti- ferromagnetic" . We also find in the latter model 
that within a semiclassical analysis a region of bistability 
arises for the first-order QPT; in the fully quantum me- 
chanical system with finite atom number, A^, signatures 
of this bistability can be identified in an atomic phase 
space distribution. We note that bistable behaviour and 
first-order non-equilibrium phase transitions have also 
been found in studies of optical bistability and resonance 
fluorescence of cooperative atomic systems fTTJ |l4l ITS] . 
However, unlike these systems, our second model involves 
direct spin-spin interaction terms and does not feature 
coherent driving of the collective atomic spin. 

A brief outline of the paper is as follows. First, in Sec. 
|TI|we briefly examine the cooperative resonance fluores- 
cence model, which exhibits a second-order QPT as the 
dissipation strength is varied, and consider the steady 
state entanglement. Then, in Sec. |III| we focus on the 
dissipative Lipkin-Meshkov-Glick (LMG) model in a pa- 
rameter regime where a second-order dissipation-driven 
QPT arises. Specifically we first present a semiclassical 
analysis of the phase transition and then consider the 
steady state entanglement behaviour across the phase 
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transition. Next in Sec. IV we present a similar anal- 
ysis for a different parameter regime where a first-order 
dissipation-driven QPT occurs in the LMG model which 
in fact exhibits bistable behaviour in the semiclassical 
steady states. Finally in Sec. [V|we will summarize our 
findings and give a brief outlook. 



II. COOPERATIVE RESONANCE 
FLUORESCENCE MODEL 

We consider here a model for cooperative resonance 
fluorescence as studied in [TTJ [121 US] , which describes a 
collection of N two-level atoms that are resonantly driven 
by a classical laser field and undergo collective sponta- 
neous emission. This system can be described by the 
following (zero-temperature) master equation, 

p = p] + ^ (2j_pj+ - j+ j_p - pj+ j_) (1) 

where Q is the strength of the coherent driving field 
and 7 is the collective spontaneous emission rate (i.e., 
7 is proportional to the atomic density) [16 . The an- 
gular momentum operators are defined in terms of the 
individual two level operators by = (1/2) cr^*\ 

J± = cr^^ , with the Pauli matrices for atomic 
spin z, and J^, = (1/2) (J+ + J_), = (-i/2)(J+ - J-). 
We note that this master equation possesses the exact 
steady state solution p/TJ [18] , 



(2) 



where J± = J± ^ iQN/{2-f). 

For a potential experimental realization of this system, 
we have in mind an ensemble of atoms coupled collec- 
tively to an optical (quantized) cavity mode and laser 
fields, which together drive Raman transitions between a 
pair of stable atomic ground states in a A-type configu- 
ration (similar to the setups described in [6^ and [T9 ). In 
particular, a pair of laser fields drive a resonant Raman 
transition between the atomic ground states to provide 
the coherent driving term in ([T]), while the cavity mode 
and another laser field drive a second, distinct Raman 
transition. In the ("bad cavity") limit where the cavity 
field decay rate is much larger than the Raman transition 
rates, the cavity mode dynamics adiabatically follows the 
atomic dynamics and can therefore be eliminated from 
the model [6 , yielding the dissipative term proportional 
to 7 in ([T]), with 7 the effective (cavity- mediated) collec- 
tive atomic spontaneous emission rate. Also note that 
in such a setup the dissipative term automatically scales 
with a factor of 1/A^ (in contrast to earlier studies [13 ), 
which allows the thermodynamic limit to be identified 
more readily and ensures that the critical point is inde- 
pendent of the system size N in this limit. 

In Sec. |II A| we first study the steady state solutions 
of the cooperative resonance fluorescence model. Then 



in Sec. II B we determine and analyze the steady state 
entanglement present in the system. 



A. Steady States 

From the above master equation we derive the fol- 
lowing semiclassical equations of motion for the compo- 
nents of the Bloch vector, X = {Jx)/j^ Y = {Jy) / ji and 
Z = {Jz)/j (see [6] and [19] for similar derivations) 



X = -fZX, 

Y = -nZ^jZY, 

Z = nY --f{X^ ^Y^), 



(3a) 
(3b) 
(3c) 



with the constraint + Y"^ + = 1 corresponding to 
conservation of angular momentum. For 7 > 7c = 1^, the 
stable steady state solutions are given by 



Zss = Xss = 0, rss = ^/7. (4) 

When 7 < 7c one finds that no (semiclassical) steady 
state solutions exists. However, the (finite-A^) master 
equation has a stable steady state solution for all 7, as 
given by Eq. ([2]), which indicates that quantum fluctua- 
tions play a crucial role in determining the state of this 
model. This was extensively discussed in earlier works 
(TTJ [12], where an effective description for A 1 was 
developed and the steady state Bloch vector components 
for 7 < 7c were shown to be 



-^ss — ^ss 



7 sin~^(7/$7) 



• (5) 



In Fig. [T] we plot the non- vanishing steady state Bloch 
vector components, as given by the above expressions, 
together with finite-A solutions (computed from numer- 
ical solution of the master equation [20 ) for comparison. 
We see that there is good agreement for sufficiently large 
values of A (in fact, the two approaches are already in 
reasonable agreement for A :^ 100). Also, in Fig. [2] we 
plot finite-A solutions for the steady state second-order 
moments, (JJ)/j^ {J^)/j^ and {J^,)/f. 
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FIG. 1: Semiclassical and asymptotic solutions (solid line), 
and finite-A steady state moments for Q = 0.2, and A = 25 
(dotted line), 50 (short dashed line), 100 (long dashed line). 
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FIG. 2: Finite- A/" steady state second-order moments for Q — 
0.2, and AT = 25 (dotted line), 50 (short dashed hne), 100 
(long dashed line). 



B. Entanglement 



for large values of N . However, in the limit of large 7 we 
run into numerical difficulties and thus we will consider 
here only the behaviour near the critical point. In the 
inset of Fig. [3] we show the behaviour of max(CR) as a 
function of N in the vicinity of the critical point (since 
at finite N the critical point depends upon N it would 
be meaningless to consider a fixed value of 7). We see 
that max(CR) continues to increase with N and appears 
to approach the asymptotic value of 1 (viz. the corre- 
sponding thermodynamic limit value, see below) with an 
approximately logarithmic scaling in N . 



0.35 




FIG. 3: Rescaled concurrence Cr for = 100 and Q = 0.2. 
Inset: Asymptotic behaviour of max(CR) as a function of N 
for — 0.2 in the vicinity of the critical point. 



We consider bipartite entanglement between individual 
atomic spins, as quantified by the rescaled concurrence, 
C'r = {N — 1)C, with C the concurrence [21 , and by the 
phase-dependent measure max{0,C(^} [22 , where 



1 



N 



(AJ; 



7V2 



(6) 



with J^p = sm{(p)Jx + cos{(p)Jy. Note that in Ref. ^ 
the rescaled concurrence was found to be related to 
through the relation Cr = max^^C^^. For our proposed 
realization of the model, the latter entanglement measure 
can in principle be determined from appropriate (quadra- 
ture variance) measurements performed on the cavity 
output field as explained in [6 . In the present model we 
find that the relation Cr = max^^ Ci^ also holds, which 
then, indirectly, enables a measurement of the rescaled 
concurrence. 

In Fig. [3] we plot the rescaled concurrence as a function 
of the dissipation strength 7 for = 100 as calculated 
numerically from the master equation (Tj). We can see 
that the entanglement peaks close to the critical point 
and then very rapidly diminishes to zero below the criti- 
cal point, in agreement with a previous study [13 . This 
behaviour can be understood by considering the entan- 
glement measure max{0,C(^} as a function of the phase 
(f and the dissipation strength 7. We find that above 
the transition, 7 > 7c, C^p is non-zero for a broad range 
of cp around cp = 7r/2. However, below the transition 
is zero for all cp and 7 as both the mean values and 
fluctuations of the components of the Bloch vector, i.e., 
{Ja)i/f and (J^)/j^ scale as f = N'^/A (see Fig.[l]and 
Fig.[2|. By using the exact steady state solution, given in 
Eq. pi) , we are able to calculate the rescaled concurrence 



In the present context it is also useful to consider a 
phase space representation of the steady state as given 
by the spin Q-function, 

QM = mv), (7) 

where 1?^) are the atomic coherent states defined by 



\ri) = 0- + \ri? 




j+rn 



b'.^)j. (8) 



with T] = e^^ tan | , where 6 and (j) correspond to spherical 
coordinates, and |j, m) are the Dicke states with m G 
[-j, -j + 1, • • • , j - 1, j] (for our system, j = N/2). 

In Fig. [4] the spin Q-function, Qs{v)^ is shown on the 
Bloch sphere for four different values of 7. We see that 
above the transition point Qs{v) is a symmetric, single- 
peaked function centered at the corresponding semiclas- 
sical amplitude. As the critical point is approached Qs{v) 
stretches and moves down to the equatorial plane. How- 
ever, once below the transition point Qs(^) remains cen- 
tered at ^ = 7r/2 for all 7 and merely continues to spread 
out in size as the fluctuations increase. Eventually, as 
7^0, Qs{v) covers the entire Bloch sphere, in agree- 
ment with the fluctuations becoming evenly distributed 
in the x, y and z directions (see Fig. [2|. 

Next, we consider the thermodynamic limit by lineariz- 
ing the quantum fluctuations around the semiclassical 
steady state solutions for A^ 1 using the Holstein- 
Primakov (HP) representation [23 . This basically in- 
volves replacing the spin ladder operators J± by bosonic 
creation (annihilation) operators (c^), where /c G {< 
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(a) 



(b) 
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FIG. 4: (Colour online) Steady state spin Q-function, Qs(^), 
on the Bloch sphere for (a) 7 = 0.05, (b) 7 — 0.15, (c) 7 = 
0.225, (d) 7 = 0.5, with AT = 50 and Q = 0.2. Note that dark 
blue corresponds to the minimum value of zero of Qs (77) while 
dark red indicates the maximum value of Qs(^)- 



III. SECOND-ORDER TRANSITION IN 
DISSIPATIVE LMG MODEL 

We now turn to the dissipative LMG model, first stud- 
ied in [6], which describes a collection of N interacting 
two-level systems in the presence of (collective) dissipa- 
tion. Specifically, as shown in [6 , by considering an en- 
semble of atoms coupled collectively to optical cavity and 
coherent laser fields, with suitably tailored Raman tran- 
sitions between a pair of atomic ground states, one may 
realize a dynamics described by the master equation 

P = -i[H^MG. p] + ^^[2 Jjp + (11) 



TV 



N 



where and r5 are tunable dissipation strengths, and 
the Hamiltonian is given by 



(12) 



with h and A tunable effective field and interaction 
strengths, respectively. 

In this section we will study this model in the regime 
A < 2/1, where a second-order phase transition occurs. In 
Sec. Ill A we analyze the steady states and non-linear 
dynamics of the semiclassical equations of motion. Then 
in Sec. IIIB we determine the steady state entanglement 
in the system. 



, >} denotes below or above the critical point, represent- 
ing the quantum fluctuations. We note that, as follows 
from the above discussion of the steady state represen- 
tation in phase space, the linearization is only possible 
for 7 > 7c, since below the critical point the state can 
no longer be described in terms of fluctuations centered 
around a semiclassical steady state solution on the Bloch 
sphere. However, above the critical point we can easily 
obtain the linearized master equation by expanding the 
angular momentum operators around the semiclassical 
steady state solution using the HP representation |23|, 
which gives 

p = 7+I)[ct]p + 7_I)[c]p+ — (2cpc-{c^p} + H.c.), 

(9) 

where c and are bosonic annihilation and creation op- 
erators, respectively, 7± = (1 — ^1 T and we have 

defined the convention D[A]p = 2ApA^ - A^Ap - pA^A. 
Using this master equation, it is straightforward to com- 
pute the rescaled concurrence analytically and we find 

Cr = 1 - ^/^^WJ^. (10) 

Away from the critical point this is in good agreement 
with the results shown for N = 100 in Fig. [Sj whilst near 
the critical point the agreement between the finite- and 
linearized results improves with increasing system size A, 
as shown in the inset of Fig. [3j 



A. Semiclassical Analysis 

1. Steady-state solutions 

The semiclassical equations of motion for the compo- 
nents of the Bloch vector, X = {Jx)/j^ Y = {Jy)/ji 
Z = {Jz)/j^ where j = A/2, derived from the above 
master equation, are given by 

X = 2hY -Ti,ZX, (13a) 
Y = -2hX ^2\ZX -T^ZY, (13b) 
Z = -2AAY' + F5(X2+r2), (13c) 

with the constraint A^ + + = 1 corresponding to 
conservation of angular momentum. 

The steady state solutions of these equations of motion 
exhibit a bifurcation at a critical dissipation strength 

VI = 2^h{\ - h), (14) 

provided X > h; also note that Fg < A. For F5 > Fg the 
stable steady-state solutions are 

Zss = -^ss = ^ss = 0, (15) 

whilst for F5 < Fg they become 

2h _ l A^-4h^ _r, 

^ss — ^ss — ±Y — 2AA — ' ~ 2/l 
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where 



the eigenvalues of M are purely real and given by 



A = A- 



4^ 



(17) 



= -r. ± r? 



(18) 



From these expressions we can see that, in the present 
regime of A < 2/i, ah the semiclassical steady state solu- 
tions vary continuously across the phase transition corre- 
sponding to a second-order phase transition. In Fig. [5] we 
illustrate this bifurcation, where, to facilitate a compar- 
ison between semiclassical and finite-A/' solutions (com- 
puted from numerical solution of the master equation), 
we plot the second-order moments (Jj) and (Jy) (since 
the finite-A^ master equation gives {Jx) = (Jy) = for 
all A), and (Jz)- We note that the two approaches are 
in reasonable qualitative agreement; we expect improved 
quantitative agreement for increasing but unfortu- 
nately we are computationally restricted from consider- 
ing much larger system sizes. This should be compared to 
the results of the model in the previous Sec. [TTj i.e., Fig.[l] 
, where the convergence between finite-A^ and asymptotic 
results was much better. Note that since A ~ a signifi- 
cant inversion of X^s or Zgg (as observed for the second- 
order transition presented in [6 ) does not occur here for 
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FIG. 5: Semiclassical (solid line) and finite-A^ steady state 
second-order moments for h = 0.2, A = 0.3, and = 0.28, 
and A^ = 25 (dotted line), 50 (short dashed line), 100 (long 
dashed line). 

Next let us briefly consider the detailed stability anal- 
ysis of the steady state solutions, which can be readily 
determined by linearizing the above non-linear equations 
of motions ( |13a[ )-( 13c) around the steady state solutions. 
The resulting linearized equations of motion can be ex- 
pressed as (X, r, Z)^ = M(X, r, Z)^ + C, where C is a 
(constant) three component vector, and we will consider 
the non-trivial eigenvalues, of the 3x3 matrix M (a 
trivial zero eigenvalue is always present due to the con- 
stant of the motion). Above the critical point, F5 > Fg, 



From this we see that the eigenvalues scale linearly with 
the dissipation (in contrast to the studies of our previous 
model [6 ) with /i+ going to zero at the critical point. 



Below the critical point, F5 
are given by 



< Fg, the eigenvalues of M 



2F5/1 



A + yX^ 



± j2(2/i2 + 2F2-AA) .(19) 



In the region < F5 < Fg, where 



(20) 



the eigenvalues are also purely real, with /i+ going to zero 
at the critical point, whilst in the region F5 < T^^ < F^ 
they become complex conjugate pairs, with an imaginary 
part that goes to zero as y^r^~^^Yl. For our character- 
istic parameters of h = 0.2 and A = 0.3 we find that 
F'^' = 0.25 and Fg = 0.28, hence F'^' < Fg. We note from 
the above expression that for the case X < 2h considered 
here, the eigenvalues vary smoothly across the critical 
point, as expected for a second-order phase transition. 



2. Time- dependent solutions 

Let us now consider numerical, time-dependent so- 
lutions o f the semiclassical equations of motion, i.e., 
of Eqs. (|13a|)-(|13c|). We have calculated the evolu- 



tion of the Bloch components X(t), Y(t), and Z{t) nu- 
merically for a uniform distribution of different initial 
states on the Bloch sphere. The resulting trajectories 
{X(t), y'(t), Z{t)} are mapped from the Bloch sphere into 
the plane using the sinusoidal projection [24 . This map- 
ping is achieved in two steps. First, the solutions for the 
Bloch components are transformed into spherical polar 
angles 0{t) and (j){t). Next, the polar angles are trans- 
formed to new coordinates U{t) and V{t) via the trans- 
formation 



uit) = m) 

vit) = e{t) - 



- tt) cos {e{t) 
7r/2. 



7r/2), 



(21) 
(22) 



In Fig. [6] we plot the trajectories using the sinusoidal 
projection for a value of F5 above the critical point (see 
the caption of Fig. 6] for numerical values of parameters). 
We can see that all initial states terminate in the unique 
steady state given by Eq. |T5| ), corresponding to the point 
/7ss = 0, = 7r/2 in Fig.pj Moreover we can see that all 
trajectories approach the stable steady state along one of 
two lines, corresponding to the stable steady state being 
a node (eigenvalues /i± being purely real). Next, we con- 
sider the case where F5 is below the critical value and plot 
the trajectories using the sinusoidal projection in Fig. [7| 
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In this case we see that different initial states terminate 
in either one of the two stable steady states given by 
Eq. (16), the corresponding values for the points Uss^ Vss 
are given in the caption of Fig. [7[ We also note that none 
of the trajectories ever terminate at Uss = 0, Fss = 7r/2, 
corresponding to the steady state (15), since this solu- 
tion is unstable in this regime. We can also see that 
the trajectories form spirals centered at the stable semi- 
classical steady states, corresponding to the eigenvalues 
being complex conjugate pairs (note that r5 < T'^ for the 
choice of parameters in Fig. IgI). 




FIG. 6: Trajectories for different initial conditions on the 
Bloch sphere using the sinusoidal projection, with h = 0.2, 
A = 0.3, and = 0.45. 
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FIG. 7: Trajectories for different initial conditions on the 
Bloch sphere using the sinusoidal projection, with h = 0.2, 
A = 0.3, and Fb = 0.15. For this choice of parameters the 
stable steady states are located at Uss — ±0.357r, Vss — 0.257r 
(as calculated from (16)). 



B. Entanglement 

We consider as before the rescaled concurrence Cr and 
the phase dependent entanglement measure max{0, 
both numerically for finite N and analytically in the ther- 
modynamic limit. Note that in this model, for finite N 
and also in the linearized analysis, we have (J^^) = 
(since there are no linear driving terms in the effective 
Hamiltonian (12) or the corresponding linearized Hamil- 
tonian [6 ), and thus C^p = 1 — {A/N){J'^). Again, we 
find that Cr 
model studied in [6j. 

In Fig. [s] (a) we plot max{0,C(^} as a function of r5 
and if ioY N = 100 with parameters corresponding to the 
second-order phase transition. We see that well above 
the transition, F5 > F^, entanglement is present for a 
broad range of angles ip. However, as the critical point 
is approached the range of angles Lp which give non-zero 
entanglement, > 0, becomes increasingly narrow. Be- 
low the transition, F5 < Fg, the region of finite C^p con- 
tinues to narrow and its maximum simultaneously shifts 
toward = as F5 decreases. 



maX(^ as for the above model and the 




FIG. 8: (a) Entanglement measure max{0,Cc^} for N — 100, 
h = 0.2, A = 0.3, and F^ = 0.28. (b) Entanglement measure 
max{0,C(^} in the thermodynamic limit, with h — 0.2, A = 
0.3, and F^ = 0.28. 

This behaviour can be explained by considering the 
spin Q-function. Figure |9] displays Qs(^) on the surface 
of the Bloch sphere for A/" = 50 and for a series of dis- 
sipation strengths F5. Above the critical point Qs{v) is 
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single-peaked and centered around the top of the Bloch 
sphere (0 = 0), with a significant rotation in a direc- 
tion between the x and y axes due to the large values 
of dissipation. Consequently we see that although C^p is 
relatively broad above the transition, stemming from the 
broad shape of the lobe, its center (where it is maximal) 
is continually shifting toward smaller values of ip as the 
rotation of Qs{v) away from the x axis decreases with 
decreasing dissipation. 

As decreases towards the critical point, Qs{v) be- 
comes increasingly elongated along a direction between 
the X and y axes, until, at the transition, it splits into 
two peaks located approximately at the two semiclas- 
sical steady state amplitudes (16). These peaks con- 



tinue to move apart in phase space as the dissipation 
strength is decreased further, approaching the corre- 
sponding dissipation- free points at 6 = 7r/2 and (j) = 
0, TT. Correspondingly, the range of (p over which C^p re- 
mains finite becomes increasingly narrow and is focussed 
around an axis perpendicular to that along which the 
two peaks lie. This narrowing of the "width" of C^p can 
be explained by noting that, since (J^^) = 0, we have 



= 1 - {4:/N){[sm{(p)Jx + cos{(p)JyY). For decreas- 
ing dissipation strength < Fg, (Jj) becomes of order 
j2 _ (^ggg p^g |5|^ gQ optimal choice of (p 

becomes more critical. 



(a) 



(b) 







(c) 




(d) 



and find that it peaks close to the critical point. In con- 
trast to the model of the previous section, significant en- 
tanglement is present on both sides of the critical point. 

We now briefiy consider the thermodynamic limit 
where we can obtain analytic results for N ^ 1 from 



the linearized HP model [6 as explained in section II B 



Since we are not directly interested in the linearized mas- 
ter equation we will refrain from quoting it here and refer 
the reader to [6]. 

In Fig. [sjb) we display C^p in the thermodynamic limit 
as a function of the dissipation strength F5 and the phase 
angle (p. We observe that above the transition, F5 > Fg, 
the behaviour is very similar to the finite- A/" result of 
Fig. (sja). However, the behaviour below the critical 
point, F5 < F^, is very different, with C^p non-zero for a 
broad range of (f due to the linearized treatment, which 
only accounts for fiuctuations around one of the two semi- 
classical steady state amplitudes (i.e., around one of the 
two lobes appearing in the spin Q-function for F5 < F^). 
Note that we can obtain plots of max{0,C(^} similar to 
Fig. [sj^a) for the region F5 < F^, but determined from 
the linearized HP model (with a finite value of A/"), by 
making a rotation back to the original coordinate system 
and then setting, by hand, (x^) = 0, to mimic an equal, 
incoherent mixture of the states associated with the two 
semiclassical amplitudes. 

Finally, in Fig.[lO]we also plot Cr as computed in the 
thermodynamic limit (solid line) from the linearized mas- 
ter equation [6 . We observe that in the linearized regime 
the peak of the concurrence is shifted below the critical 
point, whilst at the critical point there is a marked change 
in the behaviour. However, we can once again recover a 
curve within the linearized HP model that is more sim- 
ilar to the finite- A" result for F5 < F^, by following the 
procedure outlined at the end of the previous paragraph, 
which is also plotted in Fig. [To] (dot- dashed line). 

It is interesting to note that the difference in the be- 
haviour of max{0,C(^} between the finite- A^ and ther- 
modynamic limit is analogous to that observed in [6]. 
However, in [6 the rescaled concurrence peaked at the 
critical point for both the finite- A" calculations and the 
thermodynamic limit. Thus, unfortunately, the result for 
max{0, Ccp} in the thermodynamic limit does not give us 
any clues to the discrepancy observed in the rescaled con- 
currence. 



FIG. 9: (Colour online) Steady state spin Q-function, (5s(^), 
on the Bloch sphere for (a) F^ = 0.1, (b) Fb = 0.25, (c) 
Fb 0.4, and (d) Fb 1, with AT = 50, /i = 0.2, A 0.3, and 
ri — 0.28. Note that dark blue corresponds to the minimum 
value of zero of Qs{rf) while dark red indicates the maximum 
value of Qs{v)- 



In Fig. 10 we plot the rescaled concurrence calculated 
for the system sizes A" = 25, 50, 100 (set of dashed lines) 



IV. FIRST-ORDER TRANSITION IN 
DISSIPATIVE LMG MODEL 

In this section we consider the dissipative LMG model 
described by the master equation (11) in the regime 
A > 2/1, where we find that a discontinuous first-order 
transition occurs. Moreover, we find that there are in fact 
two transition points that encompass a region of bista- 
bility, the size of which increases with increasing A. In 
Sec. |III A| we analyze the steady states and non-linear 
dynamics of the semiclassical equations of motion. Then 
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FIG. 10: Rescaled concurrence Cr for = 25 (dotted line), 
= 50 (short dashed Une), N = 100 (long dashed line), 
the thermodynamic limit (solid line), and the thermodynamic 
limit in the original coordinate system (dash-dotted line, see 
text) with h = 0.2, A = 0.3, and Fl = 0.28. The vertical 
dotted line indicates the location of the critical point Fl = 
0.28. 



in Sec. |IIIB| we determine the steady state entanglement 
in the system. 



A. Semiclassical Analysis 

1. Steady state solutions 
We consider, as previously, the semiclassical equations 



of motion given by (13a)- (13c) and begin by studying 
their steady state solutions. In the region r5 > A the 



stable steady state solutions are given by (15), whilst for 

< they are given by (16). However, in the region 

< Ffo < A both steady state solutions (15) and (16) are 



in fact stable, i.e., the semiclassical system is bistable 
[25] . These two solutions each exhibit a discontinuity, 
associated with a first-order phase transition, at the crit- 
ical points and A, respectively. Note that for values 
of A not significantly larger than 2/i, the bistable region 
i s in fact very small, with F^ :^ A, and consequently 
the distinction between steady states in this region be- 
comes somewhat redundant [6]. However, in the regime 
X ^ 2h the extent of the bistable region becomes signif- 
icant and both stable steady states must be considered. 
In fact, in this situation we can expect that a complete 
description in terms of a single steady state will not be 
possible. We will focus on the regime of large A, and 
study the system primarily for the characteristic param- 
eters {h = 0.2, A = 0.75, Fa = 0.01}. 

The behaviour of the semiclassical steady state solu- 
tions as a function of F5 is illustrated in Fig. [TT] together 
with finite- A/" solutions. Outside the bistable region, con- 
vergence of the finite- solutions towards the semiclas- 
sical results is evident, while inside the bistable region 
the finite- A/" solutions appear to approach the semiclassi- 
cal branch corresponding to Eq. (16), although the rate 



of convergence with increasing A^ is clearly much slower, 
and indeed the finite- A^ curves are suggestive of some 
degree of "averaging" between the distinct semiclassi- 
cal steady states. (In fact, consideration of the spin Q- 



function later in this section will support this picture.) 

At the critical points, the semiclassical moments 
"jump" by a larger amount as A is increased; specifically 
the size of the jump at the critical points is quantified by 
A = 1 — Zss, where is given by (16). For the critical 



point F^ this is given by A = 1 — h/(X — h) and for the 
critical point A by A = 1 — 2/i/A, hence we see that at 
either critical point the largest possible jump occurs in 
the limit \ ^ h. 



0.5 



0.5 



0.7 
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FIG. 11: Stable semiclassical solutions, (16) (solid line), and 
(15]) (dash-dotted line), with h = 0.2, A = 0.75, and F^ = 
0.66. Also plotted are the finite- A/" steady state moments for 
A^ = 25 (dotted line), 50 (short dashed line), 100 (long dashed 
line). 

Now let us again consider the stability of the steady 
state solutions by examining the non-trivial eigenvalues 
li± of the matrix M describing the linearized semiclas- 
sical equations of motion. For the steady state solution 



(15), stable in the region F^ < F5, the eigenvalues are 



again given by (18), whilst for the steady state solutions 
(15), stable in the region F5 < A, they are given by (19). 



We plot these eigenvalues in Fig. [12] and note the discon- 
tinuities in (18) and (19) at F^ and A, respectively, as 
expected for a first-order transition. These discontinu- 



ities coincide with as given by (18), going to zero at 



F^, and as given by (19), going to zero at A. 



In the bistable region the behaviour of the eigenvalues 
associated with the steady states (16) is quite similar 
to that of the eigenvalues for F5 < Fg for the second- 
order transition of Sec. |III[ Note, however, that whilst 



for A < |(3+V^) (the case encountered previously in Sec. 



HI) one finds F^,' < F^, if A > |(3 



>/5) (corresponding 



12]) we find that A > F'^' 



to the case plotted in Fig 
(this was already noted in [6]) 

Finally, over essentially all of the bistable region Fg < 
F5 < A we observe that the real part(s) of the eigenvalues 
associated with the solutions (16) are smaller (i.e., more 



negative) than the largest real part of the eigenvalues 
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associated with the solutions (15). This points to the 



steady state solutions (16) being more stable than the 
solutions (15) over the majority of the bistable region, 



which is consistent with the apparent convergence of the 
finite- A/" results with increasing N towards the solutions 
(16). We note, however, that for significantly larger val- 



ues of A, where the bistable region is much larger, either 
steady state solution can be more stable dependin g on 
the choice of within the bistable region, and thus a 
convergence of the finite- results toward one of the two 
semiclassical solutions is not observed. 




FIG. 12: Eigenvalues of the linearized equations of motion, 
as given by (19) (solid line), and by (18) (dash-dotted 
line), for h = 0.2, A = 0.75, and F^ = 0.66. The inset shows 
a magnification of the real part of /i± near the critical point 
at Ffe A. 




FIG. 13: Trajectories on the Bloch sphere for different initial 
conditions, using the sinusoidal projection. Parameters are 
h = 0.2, A = 0.75, and Fb = 0.7. For this choice of parameters 
the broken phase stable steady states are located at Uss — 
zb0.567r, Vss — O.IStt (as calculated from (16)), whilst the 
normal phase steady state solution is located at Uss = 0, V^ss = 
7r/2. Note that the unstable solution described in the text is 
located at Uss — ±0.2l7r, Vss — O.SItt for the parameters 
considered. 



B. Entanglement 



2. Time- dependent solutions 

Time dependent solutions of the semiclassical equa- 
tions of motion are now examined, again using the sinu- 
soidal projection defined earlier. In particular, in Fig.p!3 
we plot trajectories {/7(t), V{t)} for a value of F5 within 
the bistable region F^ < F5 < A [26 . We see that, in this 
case, different initial states evolve to either the steady 
state (15), (corresponding to Uss = 0,^88 = 7r/2), or to 
either of the two steady states (16) (corresponding to 
Uss — ±0.567r, — O.IStt). Again we see that the way 
in which the trajectories approach the respective steady 
states is directly related to the type of eigenvalue. 

In Fig. [13] we note two distinct "gaps" centered at the 
points {U ±0.2l7r, y :^ O.SItt}. These points in fact 
correspond to unstable steady state solutions of the semi- 
classical equations of motion, given by 



2h 



A/2 _ 4/^2 



2AA' 



= ^XlZl, (23) 



where A' = A — ^/}^^^^T^. These points separate trajec- 
tories that terminate in different (stable) steady states, 
i.e., trajectories that pass "above" ("below") these points 
terminate in the steady state {Uss = 0? V^s = '7r/2} 
({/7ss = ±0.56^, ^ss^ 0.13^}). 



We now consider again the entanglement measures C^p 
and Cr, both numerically for finite N and analytically 
for ^ 1, corresponding to the linearized regime. Fig- 
ure 14 'a) shows a plot of C^n as a function of F5 and ip 



for N = 100. We see that, well above the critical value 
F5 = A, substantial entanglement is present over a broad 
range of angles (p. As F5 approaches A from above, signif- 
icant entanglement persists, but for a somewhat narrower 
range of angles ip. However, in the vicinity of F5 = A the 
entanglement diminishes rapidly for all values of (p as F5 
decreases further. 

To help understand these results we again utilize the 
atomic coherent state representation and study the spin 
Q-function. In Fig. 15 we plot Qsivj) on the Bloch sphere 
for a series of values of F5 in the vicinity of the first- 
order transition. Well above the critical point (5s(^) is 
a single-peaked function with little angular dependence. 
Correspondingly, C^p is non-zero over a broad range of 
with a maximum close io Lp = Sir /A. Note that in contrast 
to the results found in [6 , here there is a large shift of 
the optimum away from (f = 7r/2, since fiuctuations in 



both the y and x directions are significant (see Fig. 11). 

As F5 decreases towards the value A, Qsivi) becomes 
increasingly stretched along a direction between the x- 
and ?/-axes. As the critical value F5 = A is traversed, 
Qs(^) changes from a single-peaked function to a triple- 
peaked function, corresponding to the existence of three 
stable steady states in the region Fg < F5 < A. As noted 







FIG. 14: (a) Entanglement measure max{0, Cp} for N = 100, 
/i = 0.2, A = 0.75, and = 0.66. (b) Entanglement mea- 
sure max{0, Ccp} in the thermodynamic limit, with the same 
parameters, for the steady state (15). (c) Entanglement mea- 
sure max{0,C(^} in the thermodynamic limit for the steady 
state ( 16 ). 



above, the range of (f over which remains finite nar- 
rows and then drops abruptly to zero as r5 approaches 
A from above. This can be explained by noting that in 
the bistable region (J^) = 0, and thus we again have 
= 1 - {4/N){[sm{(p)Jx + cos{(p)Jy]'^). Since both 
(Jj) and {Jy) are of order = A/'^/4 in this region (see 
Fig. pTj) , this severely restricts the range of F^ and (p for 



which > 0. 



As the dissipation strength is decreased further, even- 
tually the critical point Fg is crossed and the central peak 
of Qsiv) at the top of the Bloch sphere vanishes; we then 




FIG. 15: (Colour online) Steady state spin Q-function, Qs(^), 
on the Bloch sphere for (a) F^ = 0.55, (b) F^ = 0.65, (c) F^ = 
0.705, and (d) F^ = 0.85, with iV = 50, /i = 0.2, A = 0.75 and 
F^ = 0.66. Note that dark blue corresponds to the minimum 
value of zero of Qs{ri) while dark red indicates the maximum 
value of Qs{r]). 



recover the familiar two-lobed structure associated with 
the two semiclassical steady state amplitudes. The mo- 
ments (Jj) and (Jy) are still of order = A/'^/4 in this 
region (and {J^) = 0), and consequently = for all 
(p. 

Now we briefly turn to the thermodynamic limit again 
where we can obtain analytic results for 1 from the 

linearized HP model [6^. Outside the bistable region we 
can compute the entanglement by linearizing the fluctu- 
ations about the unique stable steady state. However, 
in the bistable region we can only compute the entangle- 
ment by linearizing about one or the other of the stable 
steady states. 
In Fig. 



14 



b) we display C^p as a function of the dissi- 
pation strength F5 and the phase angle (p for the choice 
of stable steady state given by (15). We observe that, 
for F5 > F^, is non-zero over a broad range of (f cen- 
tered around cp = 37r/4, whilst near the critical point 
F^ the corresponding range of (p narrows. Outside the 
bistable region is in reasonable agreement with the 
corresponding finite- A/" results, however, as might be ex- 
pected inside the bistable region, the results differ con- 
siderably, since linearization around simply one of the 
stable steady states is inadequate. 
In Fig. 



14 



^c) we display for the choice of stable 
steady state given by ([l6|. For F5 > A, C^p is non-zero 
for a broad range of cp centered around cp = 37r/4, which 
is also in reasonable agreement with the finite- A" results. 
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However, in addition, in the bistable region < r5 < A, 
we find that is non-zero for a small range of ip. Specif- 
ically, we see that with decreasing F^, away from the 
critical point, C^p becomes broader and its center moves 
toward = tt. The existence of this second lobe and its 
behaviour can be described by considering the fluctua- 
tions in the linearized regime. In contrast to the finite- A/" 
results one finds (not actually shown) that the fluctua- 
tions below the transition are less significant, thus giving 
7^ 0. The f act that the centre of moves toward tt 
can be attributed to the dominating fluctuations in the 
X direction for F5 <C F^, as opposed to approximately 
equal fluctuations in the x and y direction for F5 :^ F^. 



Finally, in Fig. 16 we plot the rescaled concurrence 
Cr as a function of the dissipation strength F5, for the 
system sizes oi N = 25, 50, 100 (set of dashed lines) and 
in the thermodynamic limit (for each choice of stable 
steady state in the bistable region). For finite N the 
entanglement attains a peak value close to F5 = A, while 
for ^ 00 the entanglement peaks at either F5 = Fg 
or F5 = A for the different steady states, respectively. 
It is worth noting that the entanglement for the case 
corresponding to the more stable steady state, i.e., (16), 
agrees more closely with the finite- result. 




FIG. 16: Rescaled concurrence Cr for N — 2b (dotted line), 
iV = 50 (short dashed line), N = 100 (long dashed line), 
and in the thermodynamic limit [solid line for (16), and dash 
dotted line for ([l5l], with h = 0.2, A = 0.75 and F^ = 0.66. 



phase space distribution. In the dissipative LMG model 
we found that either a continuous or discontinuous phase 
transition occurs depending only on the ratio of the effec- 
tive field and interaction strengths. The steady state en- 
tanglement was analyzed in detail and the modifications 
due to strong dissipation were interpreted with the help 
of the atomic phase space distribution. In the regime of 
the first-order phase transition we showed that bistable 
behaviour can occur as evidenced by the semiclassical 
analysis; specifically we showed that whilst a linearized 
analysis is generally inadequate in this regime, one of the 
two different stable solutions tends to dominate and this 
is also reflected by finite- calculations. Finally, we have 
also briefly explained how both of these models might be 
implemented using an ensemble of atoms that interacts 
with optical cavity and laser fields, and thus how the en- 
tanglement properties might be measured via the cavity 
output field. 

In the future, it would be interesting to compare the 
phase-dependent entanglement measure considered here 
with the context-based entanglement measure studied 
in [27], where entanglement is quantified by consider- 
ing a continuous observation of the environment of a 
dissipative system. Specifically, the angle appearing in 
our phase-dependent measure can be associated with the 
phase of a local oscillator in the homodyne measurement 
considered in [27 . It would also be interesting to study 
entanglement criticality in dissipation-driven quantum 
phase transitions of other systems, such as generalized 
collective spin models [28 subjected to strong dissipa- 
tion, or in collective models with additional short-range 
interactions [29 and dissipation. Other indicative mea- 
sures of the critical behaviour might also be interesting 
to consider such as the recently proposed mixed-state 
fidelity [30 or the operator fidelity susceptibility [31]. 
Another important future direction is the simulation of 
the full quantum model for much larger system size A/", 
using, for example, quantum trajectory methods [32j, so 
that the approach towards the thermodynamic limit may 
be studied more carefully. 



V. CONCLUSION 



We have shown the presence of quantum phase tran- 
sitions occurring due to the variation of the strength of 
dissipation in two different collective spin systems. For 
the cooperative resonance fluorescence model we com- 
puted the steady state entanglement and showed how 
the results could be interpreted in terms of an atomic 
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